#ifndef SERIESAITKEN_H
#define SERIESAITKEN_H

#include "SeriesBase.h"

// A series extrapolator using the Aitken iterated delta transform

template <class T> 
class SeriesAitken : public virtual SeriesBase<T>
{
#define INITVECS(n) s.reserve(n); if (doerrcalc) { ds.reserve(n); }
public:
    SeriesAitken(uint n=20, bool _doerrcalc=true, bool _saveterms=true) : 
        SeriesBase<T>(n,_doerrcalc,_saveterms), i0(0), startsize(n) 
    { INITVECS(n); }
    SeriesAitken(const myvector<T>& v, bool _doerrcalc=true,
                 bool _saveterms=true) : 
        SeriesBase<T>(v.size(),_doerrcalc,_saveterms), i0(0),
        startsize(v.size()) 
    { INITVECS(v.size()); AddTerms(v); }
    SeriesAitken(const myvector<T>& v, const myvector<T>& verr, 
                 bool _doerrcalc=true, bool _saveterms=true) : 
        SeriesBase<T>(v.size(),_doerrcalc,_saveterms), i0(0),
        startsize(v.size()) 
    { INITVECS(v.size()); AddTerms(v,verr); }
    ~SeriesAitken() {}
#undef INITVECS

protected:
    void CalcResult(const T& term);
    const T& GetResult() const {return s[i0];}
    const myvector<T>& GetDSum() const {return ds[i0];}

private:
    myvector<T> s;  // A series of approximations to the sum
    // The latest best guess is s[i0]
    myvector<myvector<T> > ds; // ds[i][j] = d(s[i])/d(term_j)
    uint i0;
    uint startsize;
};

template <class T>
inline void SeriesAitken<T>::CalcResult(const T& term)
{
    uint n=NTerms();
    assert(s.size()+1 == n);
    if (doerrcalc) {
        assert(ds.size() == s.size());
    }

    // Update the tables for the new term
    if (n == 1) {
        s.push_back(term);
        if (doerrcalc) {
            ds.push_back(myvector<T>(1,T(1.0)));
            ds.back().reserve(startsize);
        }
        i0 = 0;
    } else {
        s.push_back(PlainSum());
        if (doerrcalc) {
            for (uint i = 0; i < n-1; i++) {
                ds[i].push_back(T(0.0));
            }
            ds.push_back(myvector<T>(n,1.));
            ds.back().reserve(startsize);
        }
        int i;
        for (i = n-3; i >= 0; i-=2) {
            T a = s[i+1]-s[i];
            T b = s[i+2]-s[i+1]-a;
            T aoverb = a/b;
            s[i] -= a*aoverb;
            if (doerrcalc) {
                for (uint j = 0; j < n; j++) {
                    T da = ds[i+1][j] - ds[i][j];
                    T db = ds[i+2][j] - ds[i+1][j] -da;
                    ds[i][j] -= 2.*da*aoverb - aoverb*aoverb*db;
                }
            }
        }
        i+=2;
        assert(i==0 || i==1);
        i0 = i;
    }
}

#endif
